EpiCovDA

Example of Model Fit and Prediction

Code

ex1 <- makeForecast(historical_window = 30, path = paste0(path,"/../"))
ex2 <- makeForecast(historical_window = 14, path = paste0(path,"/../"))

(ex1$plt + labs(subtitle = "Fitting 30 days")) + 
(ex2$plt + labs(subtitle = "Fitting 14 days")) +
  plot_layout(guides = "collect") & theme(legend.position = "bottom")

Code

baseline <- makeForecast(historical_window = 50, path = paste0(path,"/../"))
set.seed(202)
seed.list <- sample(10^5, 200)
res <- map(1:50, .f = function(j){
  fcast <- makeForecast(historical_window = 50,
                        r.seed = seed.list[j],
                        path = paste0(path,"/../"))$dat 
  data.frame(date = fcast$date,
             pred = fcast$predictedIncidence,
             iter = j)
}) |> bind_rows()

# saveRDS(list(baseline = baseline,
#              res = res),
#         "../output/AZ_06-03-2020_ens_all_hist.rds")

Example ensemble forecasts

Using most recent 14 days

Code

saved_results <- readRDS("../output/AZ_06-03-2020_ens.rds")
res           <- saved_results$res |>
  filter(!is.na(pred))
baseline      <- saved_results$baseline
mean_pred <- res |> group_by(date) |>
  summarize(pred = mean(pred))
plt <- ggplot(res) +
  geom_line(mapping = aes(x = date, y = pred, group = iter, color = "Ensemble Member"), linewidth = 1) +
  geom_line(data = mean_pred, mapping = aes(x = date, y = pred, color = "Mean Prediction"), linewidth = 1.5) +
  geom_line(baseline$dat |>
               filter(date <= baseline$forecast_date + 30), 
            mapping = aes(x = date, y = smoothedIncidence, color = "Time-Averaged Data"), linewidth = 1)+
  geom_line(baseline$dat |>
               filter(date <= baseline$forecast_date + 30), 
             mapping = aes(x = date, y = positiveIncrease, color = "Raw Data"))  +
  geom_point(baseline$dat |>
               filter(status %in% c("available")), 
             mapping = aes(x = date, y = positiveIncrease, color = "Used for Fitting")) +
  scale_color_manual(values = c("Ensemble Member" = "pink",
                                "Mean Prediction" = "deeppink2",
                                "Time-Averaged Data" = "black",
                                "Raw Data" = "black",
                                "Used for Fitting" = "deepskyblue3"))+
  theme_bw() +
  theme(legend.position = "top") +
  labs(x = "Date", y = "Incident Cases", color = element_blank()) + 
  lims(y = c(0,10000))

plt

Using 50 days of historical data

Code

saved_results_all_hist <- readRDS("../output/AZ_06-03-2020_ens_all_hist.rds")
res_all_hist          <- saved_results_all_hist$res |>
  filter(!is.na(pred))
baseline_all_hist      <- saved_results_all_hist$baseline
mean_pred_all_hist <- res_all_hist |> group_by(date) |>
  summarize(pred = mean(pred))
plt_all_hist <- ggplot(res_all_hist) +
  geom_line(mapping = aes(x = date, y = pred, group = iter, color = "Ensemble Member"), linewidth = 1) +
  geom_line(data = mean_pred_all_hist, mapping = aes(x = date, y = pred, color = "Mean Prediction"), linewidth = 1.5) +
  geom_line(baseline_all_hist$dat |>
               filter(date <= baseline$forecast_date + 30), 
            mapping = aes(x = date, y = smoothedIncidence, color = "Time-Averaged Data"), linewidth = 1)+
  geom_line(baseline_all_hist$dat |>
               filter(date <= baseline$forecast_date + 30), 
             mapping = aes(x = date, y = positiveIncrease, color = "Raw Data"))  +
  geom_point(baseline_all_hist$dat |>
               filter(status %in% c("available")), 
             mapping = aes(x = date, y = positiveIncrease, color = "Used for Fitting")) +
  scale_color_manual(values = c("Ensemble Member" = "pink",
                                "Mean Prediction" = "deeppink2",
                                "Time-Averaged Data" = "black",
                                "Raw Data" = "black",
                                "Used for Fitting" = "blue"))+
  theme_bw() +
  theme(legend.position = "top") +
  labs(x = "Date", y = "Incident Cases", color = element_blank()) + 
  lims(y = c(0,10000))

plt_all_hist

Setup for additional example

Code

new_ex <- makeForecast(state_abbr        = "CT", 
                       forecast_date     = "11-15-2020", 
                       forecast_window   = 30, 
                       historical_window = 14, 
                       r.seed            = 101,
                       path = paste0(path,"/../"))


new_ex$plt